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Today's most popular techniques for accurately calculating the dynamics of the reduced density 
operator in an open quantum system, either require, or gain great computational benefits, from 
representing the bath response function a(t) in the form a(t) = YIk p^e nRt . For some of these 
techniques, the number of terms K in the series plays the lead role in the computational cost of the 
calculation, and is therefore often a limiting factor in simulating open quantum system dynamics. 
We present an open source MATLAB program called BATHFIT 1, whose input is any spectral 
' distribution function J(oj) or bath response function, and whose output attempts to be the set of 

parameters {Pk,^k}J> =1 such that for a given value of K, the series ^JrVk^ 1 ^ i s as close as 
possible to a(t). This should allow the user to represent a(t) as accurately as possible with as few 
parameters as possible. The program executes non-linear least squares fitting, and for a very wide 
^-j. variety of forms for the spectral distribution function, competent starting values are used for these 

fits. For most forms of J(oj), these starting values, and the exact a(t) corresponding to the given 
J (to), are calculated using the recent Pade decomposition technique - therefore this program can also 
be used to merely implement the Pade decomposition for these spectral distribution functions; and 
t***~ ■ it can also be used just to efficiently and accurately calculate a(t) for any given J{oj)- The program 

C^l ' also gives the J(to) corresponding to a given a(t), which may allow one to assess the quality (in the 

a;— domain) of a representation of a(t) being used. Finally, the program can calculate the discretized 
influence functional coefficients for any J(w), and this is computed very efficiently for most forms 
of J(uj) by implementing the recent technique published in [lj]. We also provide a Mathematica 
program that can perform this last calculation, along with calculating an analytic form for these 
discretized influence coefficients, for a given analytic representation of a(t). 
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It is often useful to represent the bath response function a(t) in the form: 

K 

a (t)=Y J PRe aRt - (1) 

K 



SETTING 

The most popular open quantum system (OQS) model is currently the Feynman- Vernon model. In the Feynman- Vernon model, the 
OQS (denoted by the operator s) is coupled linearly to a set of quantum harmonic oscillators Qk'. 

H = HqQS + -ffoQS-bath + Hb at h (2) 

= Hoqs + J2 CksQ + J2 (i TO «<2- 2 + hm^lQl) ■ (3) 

K K 

In most models the quantum harmonic oscillators (QHOs) span a continuous spectrum of frequencies uj k and the strength of the 
coupling between the QHO of frequency oj and the OQS is given by the spectral distribution function J{ui): 

J H = |E^r^-^)- ( 4 ) 

K 

For the hamiltonian of the Feynman- Vernon model, the bath response function a(t) is the following integral transform of J(oj): 



1 f°° , ,{ (Puh\ 
a(t) = — J J\w) y coth y— — J cos{ujt) — i sin(wi) 
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du (5) 
e- iu]t du, J(-u) = J{oj) (6) 
dw, J(-w) = J(w), (7) 

77 J-oo 1 - exp{-pu>h) 

where, equation [7] can be written in terms of the Bose- Einstein distribution function with x = fiuih: 

,Bo S e-Emstein/ n = 1 ^_ / g v 

1 — exp(— a; J 
CONTENTS 

Non-linear least-squares fitting of the bath response function to the form YI'rPk^ 111 

The goal is to represent a(t) by the series X)i?-Pi? e ° J? * as accurately as necessary, with as low as possible a value of K. Given a 
value of K, a non-linear least-squares fitting algorithm can be used to represent a(t) by that series very accurately. MATLAB has two 
such algorithms implemented for easy use: trust region reflective, and Levenberg-Marquardt. Our MATLAB code allows the user to 
choose either method. This least-squares fitting for a chosen value of K can then be repeated for larger values of K until the resulting 
RpD has converged to satisfaction. 

Since the least-squares fitting algorithm attempts to minimize the function \a(t) — '^2%PK e ^ lRt \ with respect to the parameters 
{Pk,S1,k}, we call \a(t) — ~}2%PK eVlRt \ the objective function. 



Objective functions and starting values. 

The objective function can be obtained for any spectral distribution function from equation[S]by numerical integration, but for many 
classes of spectral distribution functions, it can be obtained much faster by evaluating analytic formulas. For the first three classes 
of spectral distribution functions below, the Pade decomposition scheme first described in 0, 0] can already represent a(t) by the 
series ^ZrVk^ 1 ^ ver y accurately with a small number of terms in the series, so one can evaluate a(t) = ^2gPg^ Kt with larger and 
larger values of K until convergence is reached, using the expressions for the parameters {p^ , fl^ } given below. For the fourth class of 
spectral distribution functions below, the objective function a(t) can be calculated very quickly with the closed form analytic expression 
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presented in that section. There are obviously other spectral distribution functions for which there are closed form expressions for 
a(t), or analytic series representations, but the forms covered in the four subsections below are enough to represent a wide range of 
physically relevant spectral distributions. When no closed form expression or analytic series is obtainable for the a(t) for a particular 
form of J(oS), one can obtain a(t) from equation [S] by numerical integration. 

Non-linear least squares fitting algorithms also typically require starting values. For the present case, these are values for the 
parameters {p R ,£l R }, such that \a(t) — YIk PR^ lRt \ i s c l° se to its global minimum with respect to {p R ,£l R }- For the first three 
spectral distribution function forms below, since the Pade decomposition scheme mentioned above already represents a(t) by the 
desired series very accurately with a small number of terms, we can use the analytic expressions for {p R , ^r} from this scheme (given 
below) as starting values for the non-linear least-squares fit. We are not aware of schemes which represent a(t) by the desired series for 
the final two forms of J(lo) listed below, so choosing good starting values for these cases will not be as simple. The final two subsections 
below discuss potentially useful strategies for guessing good starting values for these cases respectively, although these are not expected 
to bring \a(t) — '^ I R PRe Q ' itt \ as cl° se to its global minimum with respect to {p R , ^r} as are the starting values derived from the Pade 
decomposition for the first three cases. 

The expressions given below were first presented in the table 1 in the apprendix of [lj . 



Spectral distribution functions of the generalized Lorentz-Drude/Debye (gLDD) fi 
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||( 1 -E5||^8f)+ i ^ -7K-&& K€[h+l,2h] {*K,lK,"K}f =%+ . = {\ R , lR ,u-K} R=j 

Here {£^,2 js} are the \jV — 1/<sV) Pade parameters for the Bose- Einstein distribution function. Expressions for them are given in 
table H 



Spectral distribution functions of the thermally scaled generalized Lorentz-Drude/Debye (tgLDD) f 



orm 



J{w) = l - tanh V ( 2 *™ . V2 + 2 \ (11) 



K 

a(t)="£p R e Q ** (12) 

K 



where: 

PR $}r Range 



^ + if££^|^f|; -1R-&R Ke[h + l,2h] {^K^t^^i^lK^K^ 



Here {£^,2 are the \,jV — l/<yV] Pade parameters for the Fermi-Dirac distribution function. Expressions for them are given in 
table H 

Spectral distribution functions of the Meier- Tanor ( MT) form 

JM = ~ ? R + & + ^) a )k + (- - ^) 2 ) (13) 
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K 

«(*) = X>*e n ** ( 14 ) 

K 



where: 

Pk Qk Range 



i^EfeA,7/» %j^P , tfe[2/l+l,tf] {A ; .-_,}^ K {A ; .-,.J,} ; K . 

Here {£^,2 ^7} are the — \j JV\ Pade parameters for the Bose- Einstein distribution function. Expressions for them are given in 
table HI 



J(w) = Auj s e~"^ (15) 
a(t) = A5R (/T (s+1) (V s) (*(*)) + V> } (*(t) + 1)) ) + LAS* ( ff^+i ) . where > ( 16 ) 
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(i) = !(-!+ it), and, (17) 



^-)( z ) = ^Linr(z). (18) 

For s £ No (ie, if s is a non- negative integer), our definition of ip^(z) in equation [TS] is a well-known representation for the 
polygamma function, where T(z) is the well-known gamma function. For these values of s, the current versions of MATLAB and 
Mathematica have a built in implementation of the polygamma function that evaluates it in real time. However, there are various 
different generalizations of the polygamma function for negative and non-integer values of s (some popular examples can be found in 
[2] and in Mathematica's generalization of the polygamma function for s ^ No does not in fact give tp^(z) as defined in equation 
IT8l Fortunately, Paul Godfrey's implementation of the polygamma function in his MATLAB function psin(s,z) does give ip^ (z) as 
defined in equation [T5] for all s £ C with 3?(s) > 0, and the evaluation is computed in real time. This function can be found on its own 
or in Paul Godfrey's bigger special functions package both which are available for free at MATLAB Central's File Exchange. 

For this spectral distribution function, our recommendation for the starting values is less straightforward than for the above three 
cases. It may be possible to represent a(t) as a sum of complex-weighted complex exponentials analytically, but unlike the above three 
cases where these series were derived from the clever Pade decomposition, which is mathematically expected to converge with very few 
terms in the series, we do not know of any such series that represents a(t) for this J(oj) such that convergence will be achieved with 
very few exponentials. Therefore, we recommend that the objective function a(t) is calculated using the analytic formula presented 
above, and that the starting parameters are chosen based on observing some of its properties, such as its frequency of oscillation and 
damping rate. Guidelines for choosing starting parameters this way are presented in the subsection below. 

Other spectral distribution function forms When J{u>) is only known in a form such that a(t) is not easily represented by equation 
[1] (as for the case in subsection,, or if J(oS) is only known numerically), it can be much harder to find good starting values for the 
non-linear least squares fit. It might be useful to fit J(u>) to a form for which good starting parameters are easy to choose, such as 
the first three forms presented above, but if this is not easy, one can calculate the objective function a(t) by numerical integration of 
equation [SJ and as mentioned towards the end of subsection „ one can use properties of a(t) as a guide to choose starting values for 

Our recommendation for the starting values for {p^} is based on the fact that at t = 0, equation [T] gives us the relation: 

K 

a{0)=Y,PK (19) 

K 

We recommend to first attempt to fit one term {p\e lf ) to the objective function a(t), withpx = a (0) according to equation [19] (bear 
in mind that will then be because sin(o; ■ 0) = will nullify the right side of equation [5]) . The strating parameter for 51i can 

then be chosen by looking at a plot of the objective function a(t) and comparing it to the expressions: 
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R(a(t)) = Pl e R{ni) cos(SJ(fii)t) , (20) 

3r(a(i)) = pie S(0l) sin(3(f2i)t) . (21) 

Based on equation[S]we see that as the temperature gets higher, 3i(a(i)) becomes more and more different from 9(a(t)), and therefore 
more than one damping rate (5R(f2^-)) and more than one angular frequency (3(f2^-)) will be needed for a good fit. However, as a crude 
estimate, we can choose the starting value of 3(f2i) to be an average of the angular frequencies (f2 = ^r, T =period of oscillations) 
of the real and imaginary parts of the objective function a(t). Since the effect of J(w) on a(t) will cause the frequency of oscillations 
of a(t) to be less and less sinusoidal over time, we recommend to estimate T of each complex component of a(t) based on the time it 
takes that component to get to the first quarter (or half) of its first oscillation. 

With this starting value for 3(f2i), the starting value of SR(f2i) can then be chosen to be an average of the damping rates of the real 
and imaginary parts of the objective function a(t). For the real part of a(t), the damping rate JJ(Qi,sr) can be estimated by observing 
the time i*- at which Sft(a(i)) attains its first minimum osr, and then solving the equation: 



ajf = exp ( 3?(Qi sft)t -a- ) cos(9 ! (f2i)i * ) , whose solution is (22) 



In 



pi cos(G(f2i)t 5 ) 

R(fi liSl ) = ^ (23) 

t -4" 



which is an estimate of ^(fii,^ ). 

Similarly, the damping rate of the imaginary part of a(t), which we denote by ^(Oi^), can be estimated bj0: 



In 



pi sin(3(fii)t ) 



= — ; — • (24) 

t 

5R(fii) is then estimated as an average of K(J7i,sr) and ^(fii^). 

We can then fit the one term function p\e l( to the objective function, and then use the resulting fitted values of p\ and Sli as 
starting values for a fit of the two term series p\e l4 + P2e to the objective function. Since P2 and O2 are expected to be of the same 
order of magnitude as piand Ox respectively, we can multiply piand f2i by random numbers near 1 in order to get crude estimates of 
suitable starting values for £»2and O2 respectively. In our MATLAB program, these random numbers are obtained using (1+RANDN) 
in MATLAB, rather than RAND or RANDN, so that these random numbers are more likely to be close to 1 rather than to 0. The 
resulting values of pi,f2i,p2 and O2 from this fit can then be used as starting values for a three term fit, with starting values for p% 
and O3 chosen as pi and O2 were for the two term fit. 

Fits to series with larger numbers of terms can be done in a similar way, although when there are many terms in the series, it may 
be more appropriate (based on equation [T9|) for the starting value of p\ to be a (°)/K instead of just a(0), though the other starting 
values would likely also have to be adjusted, which would not be straightforward. 

Note about implementation: Constraints, scaling & step sizes, and weights 

Putting constraints on the fitting parameters can help speed up the fits, and can prevent the fitting program from getting lost in 
a far from optimal local minimum, or the fitted values form becoming grossly unphysical. For example, it helps to implement the 
constraint 5R(0^ ) < 0, so that a(t) is not likely to diverge. Adding this as a constraint to the fit will prevent the fitting program from 
bothering to try values that we know are expected to give bad results, and will therefore speed up the fitting caluclation, and could 
potentially prevent the fitting program from getting 'lost' in a region far from the desired global minimum. We have implemented these 
constraints in BATHFIT's fitting routine. When K is small, it may also help to implement the constraint p\ > 0, since according to 
equation [T9l and l5| if K = 1, pi > 0. 

Since MATLAB's fitting routines are most easily implemented when the step sizes are the same size, we scale the objective function's 
range to be between and 1, and its domain is mapped to t € [0, 1], which significantly helps preventing the fitting program from 
getting lost in regions far from the desired global minimum. After the fit is complete, we then scale {pk^k} back to SI units. 

We also provide the user with the option of assigning weights to each datapoint of the objective function. For example, if the user 
requires the first lOOfs of a(t) to be represented very accurately, and does not require a(t) to be represented very accurately after 500fs, 
the user can assign weights accordingly. Often this can also help to get a graphically better fit, since the local minimum found by 
BATHFIT won't always be the closest to the global minimum, or even if the global minimum is attained, it may have been influenced 
too strongly by certain parts of a(t) which are not so important. 



1 Since this expression uses the amplitude pi which is derived from 5R(o(t)), and !R(o(t)) is expected to deviate more and more from 3(a(t)) as the temperature 
is increased (based on equation [5j, this estimate is expected to be best at low temperatures. 
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Analytic coefficients for the discretized influence functional 

Once a(t) is represented in the form of equation [TJ the DIFs (discretized influence functional coefficients, which are required if using 
a Feynman integral to calculate the RpD of open quantum systems modeled by the Feynman- Vernon model), can be calculated very 
quickly by the analytic formulas which were presented in These formulas are listed again below, along with analogous formulas 
for spectral distribution functions of the form J(oj) = Aoj s e~^' Uc ' . For this latter form for the spectral distribution function, a(t) can 
once again be fitted to the form of equation [TJ and the DIFs can therefore be calculated with the first set of formulas below, but since 
we have an exact analytic form for a(t) that can easily be integrated with respect to t, we have also presented analytic forms for the 
DIFs for this form of J(lo), which do not require a(t) to first be fitted to the form of equation [JJ 

Spectral distribution functions with bath response functions of the form a(t) = Y^,k PR eClRt 



Trotter splitting 



Strang splitting 
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4 E §T sinh 2 (r^At/ 2 ) e ^(fc-fc')A* ; o < fc' < fc < TV , and 



2 §jt ( s™H Q K At Me"** t/2 - ^0, R At ) , < k < N . 



-i K 



K=l 



Vno 



Voo = Vnn 



4 ^g Le ^(*--/ 2 ) (sinh 2 (^A t/4 )) , 

K 

2 E nT (e lRAt/i sinh("*At/ 4 ) _ Atn^ , 



= 1 K 



K 



4^ s inh(^KAt/ 2 ) sinh (fi*At/ 4 ) e ^(fcA t -Ay 4) 

j=l K 



K 



VNk = 



sinh(«*At/ 2 ) S inh (n* At/ 4 ) e n*(t-fcAt-Ay 4 ) 

j=l K 



(25) 
(26) 

(27) 
(28) 
(29) 
(30) 



Spectral distribution functions with exponential cut-offs 
QUAPI 

When the bath is nearly adiabatic, it is helpful to rewrite equation [3] as d, 



H — HoQS — ^displacement + HoQS — bath ~h ^bath ^displacement 
2 2 

= ffoQS - E + E c ^ + E + hm^lQl) + E 

ft K ft K/ 

= HqQS, displaced + E C * S Q + -^bath, displaced- 



2m K uj% 



(31) 
(32) 

(33) 



^displacement is called the "counter term", and can also be represented in terms of the spectral distribution function by recognizing that 
when the QHOs span a continuous spectrum of frequencies lo k , we have the relation (remembering equation [4]) : 



E 



1 f°° J(u) 



duj . 



(34) 



A Feynman integral used to calcualte the RpD for an OQS denoted by the hamiltonian Hoqs, displaced, is called a QUAPI, which 
stands for Quasi- Adiabatic Propagator FeynmarJl Integral. For QUAPI calculations, all r\ coefficients remain the same as for a 



2 The term 'path integral' is used more commonly than 'Feynamn integral' here, but this term is ambiguous. Currently, the first result on the search engine 
at www.google.com, when the search query 'path integral' is entered, is a Wikipedia page that currently links to three different meanings of the word 'path 
integral': (f ) line integral, (2) functional integration, and (3) path integral formulation. Only the third of these is unambiguously the Feynman integral 
discussed in this paper. The 'line integral' is an integral over a path, rather than over a set of paths; and the term 'functional integral' can refer to at least 
three types of functional integrals: (f) the Wiener integral, (2) the Levy integral, and (3) the Feynman integral. 
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feynman integral used to calculate the RpD for an OQS whose hamiltonian has not been displaced according to equations [c 
except for the r) kk t coefficients when k = k' (Vfc £ [0, N]), which change to (see supplementary material of [f|): 

% TAPI = mk ^fMd,, for fc £ [0, N] (35) 
Tvk Jo uj 

= ** + ^A. (36) 

where in the last line we have defined the "bath reorganization energy" by A. An analytic expression exists for A for most forms of 
J(oj) presented in this paper: 

JH t 

TV/fX 1 7TOJ Aft 7T 2 V^'l \ 

2 ^ h (72 + (-+- h ) a )(72 + ("-^) a ) s^Fh^J^^ 



APPENDIX 



Table I. [,yK — l,,yV] Pade parameters for the Bose-Einstein and Fermi-Dirac distribution functions (first presented in 0] and in more detail in 
Q). {Ej;,£j?} can be calculated easily for all jY , for arbitrary values of jY in our open source MATLAB program that supplements this paper. 
The matrix A is a 2jV X2^Y , and A is a 2,yV -\y,2,jV -\. The indices ,JV run from 1 to half the number of non-zero eigenvalues of the corresponding 
matrix. 

(£ | . — £ , ); and for matrices with an odd number of eigenvalues, the unpaired eigenvalue will always be 0. 

/3ft-eigcnvalues( A) fill- eigenvalues (A) 

Bose-Einstein Fermi-Dirac 

A - S m,n±l A 
1 ^mn — / ■ J *-mn — 



^/(2m + l)(2n + l) m " ^/(2m-l)(2n-l) 

A - *m,n±l A _ <*m,»±l 

y(2m + 3)(2n+3) (2m+l) (2n + l) 

-./K — \yf + 2 ,yY I Tv" <h -e 2 ) "-^ — + 2^ JTT* (e 2 -e 2 i 
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